#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
3×3 采样点 → 计算 f、fx、fy、fxy → 与 AINV 相乘得到双立方系数
并绘制 5×5 细分网格的插值热力图
"""

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt

# -------------------------------------------------
# 1️⃣ 示例 3×3 数据（请自行替换为真实像素值）
f_grid = np.array([
    [10., 12., 14.],
    [16., 18., 20.],
    [22., 24., 26.]
], dtype=float)          # shape (3,3)

# -------------------------------------------------
# 2️⃣ 中心差分求偏导
fx = (f_grid[1, 2] - f_grid[1, 0]) / 2.0          # ∂f/∂x at (1,1)
fy = (f_grid[2, 1] - f_grid[0, 1]) / 2.0          # ∂f/∂y at (1,1)
fxy = (f_grid[2, 2] - f_grid[2, 0] - f_grid[0, 2] + f_grid[0, 0]) / 4.0  # ∂²f/∂x∂y
f_center = f_grid[1, 1]

# -------------------------------------------------
# 3️⃣ 把 3×3 填入 4×4（左下角对齐），缺失的第 4 列/行用最右/上边界复制
F = np.zeros((4, 4), dtype=float)
F[:3, :3] = f_grid
F[3, :3] = f_grid[2, :]          # 第 4 行复制最上面一行
F[:3, 3] = f_grid[:, 2]          # 第 4 列复制最右边一列
F[3, 3] = f_grid[2, 2]           # 右下角复制右上角

FX = np.zeros((4, 4), dtype=float)
FY = np.zeros((4, 4), dtype=float)
FX[1, 1] = fx
FY[1, 1] = fy

# -------------------------------------------------
# 4️⃣ 按 AINV 列顺序展开为 16 维向量 c
c = np.hstack([
    F[0, 0], F[1, 0], F[2, 0], F[3, 0],   # f00, f10, f20, f30
    FX[0, 0], FX[1, 0], FX[2, 0], FX[3, 0],   # fx00 … fx30
    F[0, 1], F[1, 1], F[2, 1], F[3, 1],   # f01, f11, f21, f31
    FY[0, 1], FY[1, 1], FY[2, 1], fxy      # fy01, fy11, fy21, fxy11
])

# -------------------------------------------------
# 5️⃣ 逆矩阵 AINV（你提供的 16×16）
AINV = np.array([
    [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [-3, 3, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [2, -2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2, -1, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0],
    [-3, 0, 3, 0, 0, 0, 0, 0, -2, 0, -1, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, -3, 0, 3, 0, 0, 0, 0, 0, -2, 0, -1, 0],
    [9, -9, -9, 9, 6, 3, -6, -3, 6, -6, 3, -3, 4, 2, 2, 1]
], dtype=float)

# -------------------------------------------------
# 6️⃣ 计算双立方系数向量 a = AINV @ c
a = AINV @ c
A_coeff = a.reshape(4, 4)   # 4×4 系数矩阵

print("=== 双立方系数矩阵 (4×4) ===")
print(A_coeff)

# -------------------------------------------------
# 7️⃣ 细分网格（5×5）并计算插值
def bicubic_interp(u, v, coeff):
    """u、v 为局部坐标 (0~1)，coeff 为 4×4 系数矩阵"""
    U = np.array([1, u, u**2, u**3])
    V = np.array([1, v, v**2, v**3])
    return U @ coeff @ V.T

grid_u = np.linspace(0, 1, 5)
grid_v = np.linspace(0, 1, 5)
interp_vals = np.zeros((5, 5))
for i, u in enumerate(grid_u):
    for j, v in enumerate(grid_v):
        interp_vals[i, j] = bicubic_interp(u, v, A_coeff)

print("\n=== 5×5 细分网格插值值 ===")
print(interp_vals)

# -------------------------------------------------
# 8️⃣ 可视化（热力图）
plt.figure(figsize=(6, 5))
cmap = plt.get_cmap('viridis')
im = plt.imshow(interp_vals, extent=[0, 1, 0, 1],
                origin='lower', cmap=cmap, interpolation='nearest')
plt.title('Bicubic interpolation of 3×3 grid (5×5 sample)')
plt.xlabel('u (normalized x)')
plt.ylabel('v (normalized y)')
plt.colorbar(im, label='Interpolated value')
plt.grid(False)

# 保存本地文件（可自行打开查看）
plt.savefig('bicubic_interpolation.png', dpi=150, bbox_inches='tight')
plt.close()

